#pragma once

#define lambdaQCD 0.250
#define Z3 1.20206

double q_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2);
double S_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2);
double P_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2);
double T_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2);
void Z4F_evolution_to_RGI(O55 wm1,int Nf,double ainv,double mu_2);
double q_evolution_to_RGI(int Nf,double ainv,double mu_2);
double alphas(int Nf,double ainv,double mu);
void RI_evolution(O55 wm1,int Nf,double alm);

void evolve_Z4F_to_RGI( O55 ZpPC_RGI, O55 ZpPC, double P2, int Nf, double ainv)
{
	// Mauro PhD thesis
	// A.36 (O^RGI)^T=O^T(mu) w^-1(mu)
	// w^-1(mu) from A.35 
	O55 wm1;
	O55 ZpPC_T;
	O55 ZpPC_RGI_T;
	Z4F_evolution_to_RGI(wm1, Nf,ainv,P2);

	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC_RGI_T[iop1][iop2]=0;
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC_T[iop1][iop2]=ZpPC[iop2][iop1];
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) {
		for (int iop=0; iop<5; iop++) ZpPC_RGI_T[iop1][iop2]= ZpPC_RGI_T[iop1][iop2]+ZpPC_T[iop1][iop]*wm1[iop][iop2];}
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC_RGI[iop1][iop2]=ZpPC_RGI_T[iop2][iop1];
}
double evolve_Zq_to_RGI (double Zq,  double P2, int Nf, double ainv)
{
	double cq;
	cq=q_evolution_to_RGI(Nf,ainv,P2);
	double Zq_RGI;	
	Zq_RGI=Zq/cq;
	return Zq_RGI;
}
void evolve_to_RIp_ainv(double *Zq_RIp_ainv, double *ZO_RIp_ainv, double Zq,double *ZO, double P2, int Nf, double ainv)
{

  // scale evolution coefficients
  // from mu^2=p^2 to mu=1/a
  double cq;
  double *cO=malloc(sizeof(double)*(5));

  // Note that ZV  ZA are RGI because they're protected by the WIs
  cq=q_evolution_to_RIp_ainv(Nf,ainv,P2);
  cO[0]=S_evolution_to_RIp_ainv(Nf,ainv,P2); //S
  cO[1]=1.; //V
  cO[2]=P_evolution_to_RIp_ainv(Nf,ainv,P2); //P
  cO[3]=T_evolution_to_RIp_ainv(Nf,ainv,P2); //T
  cO[4]=1.; //A

  Zq_RIp_ainv[0]=Zq/cq;
  for (int igamma=0; igamma<5; igamma++) ZO_RIp_ainv[igamma]=ZO[igamma]/cO[igamma];
}

///////////////////////////////////
// evolution from mu=p to mu0=1/a
// Z(mu0)=Z(mu) c(mu0)/c(mu)
// def: c=c(mu)/c(mu0)
// -> Z(mu=1/a) = Z(mu) /c
//////////////////////////////////
double q_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2)
{
	double cmu, cmu0; // c=cmu/cmu0
	//mu_2=a2*p2 (adimensional quantity)
	//mu0_2=a2*(1/a^2)=1
	double mu0_2=1; 	
		
	// alphas @ NNLO
	double alm, al0;
	alm=alphas(Nf,ainv,mu_2)/(4*PI);
	al0=alphas(Nf,ainv,mu0_2)/(4*PI);

	////////////////////////////////
        // N3LO FORMULA
	// Assuming landau gauge
	///////////////////////////////////
       if(Nf==2){
        cmu = 1. + 2.03448 * alm + 35.9579 * pow(alm,2) + 1199.16 * pow(alm,3);
        cmu0 = 1. + 2.03448 * al0 + 35.9579 * pow(al0,2) + 1199.16 * pow(al0,3);
       }if(Nf==0){
        cmu = 1. + 2.0303 * alm + 42.1268 * pow(alm,2) + 1728.43 * pow(alm,3);
        cmu0 = 1. + 2.0303 * al0 + 42.1268 * pow(al0,2) + 1728.43 * pow(al0,3);
       }if(Nf==4){
        cmu = 1. + 2.4000 * alm + 29.6724 * pow(alm,2) + 719.141 * pow(alm,3);
        cmu0 = 1. + 2.4000 * al0 + 29.6724 * pow(al0,2) + 719.141 * pow(al0,3);
       }



	return cmu/cmu0;
}
double S_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2)
{
        double cmu, cmu0; // c=cmu/cmu0
        //mu_2=a2*p2 (adimensional quantity)
        //mu0_2=a2*(1/a^2)=1
        double mu0_2=1;
        
        // alphas @ NNLO
        double alm, al0;
        alm=alphas(Nf,ainv,mu_2)/(4*PI);
        al0=alphas(Nf,ainv,mu0_2)/(4*PI);

        ////////////////////////////////
        // N3LO FORMULA
        // Assuming landau gauge
        ///////////////////////////////////

      if(Nf==2){
        cmu = pow(alm,-12./29) * (1. - 8.55727 * alm - 125.423 * pow(alm,2) -
             3797.71 * pow(alm,3));

        cmu0 = pow(al0,-12./29) * (1. - 8.55727 * al0 - 125.423 * pow(al0,2) -
             3797.71 * pow(al0,3));
      }if(Nf==0){
        cmu = pow(alm,-4./11) * (1. - 8.08264 * alm - 151.012 * pow(alm,2) -
             5247.93 * pow(alm,3));

        cmu0 = pow(al0,-4./11) * (1. - 8.08264 * al0 - 151.012 * pow(al0,2) -
             5247.93 * pow(al0,3));
      }if(Nf==4){
        cmu = pow(alm,-12./25) * (1. - 9.38987 * alm - 96.2883 * pow(alm,2) -
             2403.82 * pow(alm,3));

        cmu0 = pow(al0,-12./25) * (1. - 9.38987 * al0 - 96.2883 * pow(al0,2) -
             2403.82 * pow(al0,3));
      }



        return cmu/cmu0;
}

double P_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2)
{
        double cmu, cmu0; // c=cmu/cmu0
        //mu_2=a2*p2 (adimensional quantity)
        //mu0_2=a2*(1/a^2)=1
        double mu0_2=1;
        
        // alphas @ NNLO
        double alm, al0;
        alm=alphas(Nf,ainv,mu_2)/(4*PI);
        al0=alphas(Nf,ainv,mu0_2)/(4*PI);

        ////////////////////////////////
        // N3LO FORMULA
        // Assuming landau gauge
        ///////////////////////////////////
      if(Nf==2){
        cmu = pow(alm,-12./29) * (1. - 8.55727 * alm - 125.423 * pow(alm,2) -
             3797.71 * pow(alm,3));

        cmu0 = pow(al0,-12./29) * (1. - 8.55727 * al0 - 125.423 * pow(al0,2) -
             3797.71 * pow(al0,3));
      }if(Nf==0){
        cmu = pow(alm,-4./11) * (1. - 8.08264 * alm - 151.012 * pow(alm,2) -
             5247.93 * pow(alm,3));

        cmu0 = pow(al0,-4./11) * (1. - 8.08264 * al0 - 151.012 * pow(al0,2) -
             5247.93 * pow(al0,3));
      }if(Nf==4){
        cmu = pow(alm,-12./25) * (1. - 9.38987 * alm - 96.2883 * pow(alm,2) -
             2403.82 * pow(alm,3));

        cmu0 = pow(al0,-12./25) * (1. - 9.38987 * al0 - 96.2883 * pow(al0,2) -
             2403.82 * pow(al0,3));
      }




        return cmu/cmu0;
}
double T_evolution_to_RIp_ainv(int Nf,double ainv,double mu_2)
{
        double cmu, cmu0; // c=cmu/cmu0
        //mu_2=a2*p2 (adimensional quantity)
        //mu0_2=a2*(1/a^2)=1
        double mu0_2=1;
        
        // alphas @ NNLO
        double alm, al0;
        alm=alphas(Nf,ainv,mu_2)/(4*PI);
        al0=alphas(Nf,ainv,mu0_2)/(4*PI);

        ////////////////////////////////
        // N2LO FORMULA
        // Assuming landau gauge
        ///////////////////////////////////

       if(Nf==2){
        cmu = pow(alm,4./29) * (1. + 2.66852 * alm + 47.9701 * pow(alm,2));

        cmu0 = pow(al0,4./29) * (1. + 2.66852 * al0 + 47.9701 * pow(al0,2));
       }if(Nf==0){
        cmu = pow(alm,4./33) * (1. + 2.53260 * alm + 57.8740 * pow(alm,2));

        cmu0 = pow(al0,4./33) * (1. + 2.53260 * al0 + 57.8740 * pow(al0,2));
       }if(Nf==4){
        cmu = pow(alm,4./25) * (1. + 2.91662 * alm + 37.9471 * pow(alm,2));

        cmu0 = pow(al0,4./25) * (1. + 2.91662 * al0 + 37.9471 * pow(al0,2));
       }



        return cmu/cmu0;
}

// evolve from mu=p  to RGI
void Z4F_evolution_to_RGI(O55 wm1,int Nf,double ainv,double mu_2)
{
	double alm;
	alm=alphas(Nf,ainv,mu_2);
	RI_evolution(wm1, Nf,alm);
}
void RI_evolution(O55 wm1,int Nf,double alm)
{
double beta_0 = (11.*Nc-2.*Nf)/3.;
wm1[0][0]=(4 + (alm*((-17397 + 2*(1035 - 52*Nf)*Nf)/(6.*pow(33 - 2*Nf,2)) + log(256)))/PI)/(4.*pow(alm,2/beta_0));
wm1[1][0]=0;
wm1[2][0]=0;
wm1[3][0]=0;
wm1[4][0]=0;
wm1[0][1]=0;
wm1[1][1]=(96 + (alm*(pow(alm,9/beta_0)*pow(33 - 2*Nf,2)*(-15*(13 + 32*log(2)) + 8*Nf*(1 + log(4))) -  2*(-270*(2813 + 1936*log(2)) -  6*pow(Nf,2)*(1193 + 496*log(2) + 88*log(4)) +  9*Nf*(13073 + 8008*log(2) + 484*log(4)) +    8*pow(Nf,3)*(17 + log(256)))))/(pow(33 - 2*Nf,2)*(-30 + Nf)*PI))/(96.*pow(alm,1/beta_0));
wm1[2][1]=(-96*(-1 + pow(alm,9/beta_0)) +  (alm*(-8100*(-115 + 1936*log(2)) -  54*pow(Nf,2)*(10353 + 17200*log(2)) +        270*Nf*(12843 + 28336*log(2)) +  pow(alm,9/beta_0)*(-3 + Nf)*(pow(Nf,2)*(89364 - 4032*log(2)) +     135*(90203 - 3872*log(2)) + 36*Nf*(-44915 + 2244*log(2)) +  32*pow(Nf,3)*(-51 + log(4))) +  12*pow(Nf,3)*(2549 + 3504*log(2) + 8*log(4)) -  16*pow(Nf,4)*(33 + 10*log(16))))/ (pow(33 - 2*Nf,2)*(-30 + Nf)*(-3 + Nf)*PI))/(144.*pow(alm,1/beta_0));
wm1[3][1]=0;
wm1[4][1]=0;
wm1[0][2]=0;
wm1[1][2]=(pow(alm,(8 + beta_0)/beta_0)*(195 + 480*log(2) - 8*Nf*(1 + log(4))))/ (64.*(-30 + Nf)*PI);
wm1[2][2]=(pow(alm,8/beta_0)*(96 +        (alm*(135*(-90203 + 3872*log(2)) +           4*Nf*(404235 - 20196*log(2) + 3*Nf*(-7447 + 336*log(2)) - 8*pow(Nf,2)*(-51 + log(4)))))/(pow(33 - 2*Nf,2)*(-30 + Nf)*PI)))/96.;
wm1[3][2]=0;
wm1[4][2]=0;
wm1[0][3]=0;
wm1[1][3]=0;
wm1[2][3]=0;
wm1[3][3]=(24*(-16 + sqrt(241) +  (16 + sqrt(241))*pow(alm,(2*sqrt(241))/(3.*beta_0))) +  (alm*(6*(-16 + sqrt(241) +   (16 + sqrt(241))*pow(alm,(2*sqrt(241))/(3.*beta_0)))*((1740825 + 8*Nf*(-267153 +  2*Nf*(17222 + Nf*(-1020 + 19*Nf))))/(2.*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)) + (10*log(2))/3.) - ((-1 + pow(alm,(2*sqrt(241))/(3.*beta_0)))*(-3727125 - 6552504*Nf + 1032800*pow(Nf,2) - 63168*pow(Nf,3) + 1168*pow(Nf,4) + 28*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*log(2)))/(pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf))))/PI)/(48.*sqrt(241)*pow(alm,(1 + sqrt(241))/(3.*beta_0)));
wm1[4][3]=(-72*(-1 + pow(alm,(2*sqrt(241))/(3.*beta_0))) + (alm*((-16 + sqrt(241) +  (16 + sqrt(241))*pow(alm,(2*sqrt(241))/(3.*beta_0)))*(-300015 + 370872*Nf - 39328*pow(Nf,2) + 3648*pow(Nf,3) -             80*pow(Nf,4) + 28*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*log(2)) -   (-1 + pow(alm,(2*sqrt(241))/(3.*beta_0)))*(-17978175 + 22206168*Nf - 3143648*pow(Nf,2) + 194496*pow(Nf,3) - 3664*pow(Nf,4) + 116*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*log(2))))/(pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*PI))/(144.*sqrt(241)*pow(alm,(1 + sqrt(241))/(3.*beta_0)));
wm1[0][4]=0;
wm1[1][4]=0;
wm1[2][4]=0;
wm1[3][4]=(360*(-1 + pow(alm,(2*sqrt(241))/(3.*beta_0))) +     (alm*(-6*(-16 + sqrt(241))*(16 + sqrt(241))*(-1 + pow(alm,(2*sqrt(241))/(3.*beta_0)))*((1740825 + 8*Nf*(-267153 + 2*Nf*(17222 + Nf*(-1020 + 19*Nf))))/ (2.*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)) +     (10*log(2))/3.) +     ((16 + sqrt(241) + (-16 + sqrt(241))*pow(alm,(2*sqrt(241))/(3.*beta_0)))*(-3727125 - 6552504*Nf + 1032800*pow(Nf,2) -   63168*pow(Nf,3) + 1168*pow(Nf,4) +     28*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*log(2)))/(pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf))))/PI)/(48.*sqrt(241)*pow(alm,(1 + sqrt(241))/(3.*beta_0)));
wm1[4][4]=(72*(16 + sqrt(241) +  (-16 + sqrt(241))*pow(alm,(2*sqrt(241))/(3.*beta_0))) +    (alm*(-((-16 + sqrt(241))*(16 + sqrt(241))*(-1 + pow(alm,(2*sqrt(241))/(3.*beta_0)))*(-300015 + 370872*Nf - 39328*pow(Nf,2) + 3648*pow(Nf,3) -              80*pow(Nf,4) +     28*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*log(2))) +   (16 + sqrt(241) + (-16 + sqrt(241))*pow(alm,(2*sqrt(241))/(3.*beta_0)))*          (-17978175 + 22206168*Nf - 3143648*pow(Nf,2) +           194496*pow(Nf,3) - 3664*pow(Nf,4) + 116*pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*log(2))))/(pow(33 - 2*Nf,2)*(125 + 4*(-33 + Nf)*Nf)*PI))/   (144.*sqrt(241)*pow(alm,(1 + sqrt(241))/(3.*beta_0)));
}

double alphas(int Nf,double ainv,double mu2)
{

   int CF;
   CF = (Nc*Nc-1.)/(2.*Nc);

   double lam0,L2,LL2,b1,b2,b3;
   double als0, als1, als2, als3;
   double beta_0,beta_1,beta_2,beta_3;

      beta_0 = (11.*Nc-2.*Nf)/3.;

      beta_1 = 34./3.*pow(Nc,2) - 10./3.*Nc*Nf-2*CF*Nf;
      beta_2 = (2857./54.)*pow(Nc,3) + pow(CF,2)*Nf -
             205./18.*CF*Nc*Nf -1415./54.*pow(Nc,2)*Nf +
             11./9.*CF*pow(Nf,2) + 79./54.*Nc*pow(Nf,2);

      beta_3 = (150653./486. - 44./9.*Z3)*pow(Nc,4) +
            (-39143./162. + 68./3.*Z3)*pow(Nc,3)*Nf +
            (7073./486. - 328./9.*Z3)*CF*pow(Nc,2)*Nf +
            (-2102./27. + 176./9.*Z3)*pow(CF,2)*Nc*Nf +
             23.*pow(CF,3)*Nf + (3965./162. + 56./9.*Z3)*pow(Nc,2)*pow(Nf,2) +
            (338./27. - 176./9.*Z3)*pow(CF,2)*pow(Nf,2) +
            (4288./243. + 112./9.*Z3)*CF*Nc*pow(Nf,2) + 53./243.*Nc*pow(Nf,3) +
             154./243.*CF*pow(Nf,3) +
            (-10./27. + 88./9.*Z3)*pow(Nc,2)*(pow(Nc,2)+36.) +
            (32./27. - 104./9.*Z3)*Nc*(pow(Nc,2)+6)*Nf +
            (-22./27. + 16./9.*Z3)*(pow(Nc,4) - 6.*pow(Nc,2) + 18.)/pow(Nc,2)*pow(Nf,2);

      b1=beta_1/beta_0/4./PI;
      b2=beta_2/beta_0/16./pow(PI,2);
      b3=beta_3/beta_0/64./pow(PI,3);

      lam0=lambdaQCD/ainv;

      L2   = log( mu2/(pow(lam0,2) ) );
      LL2  = log( L2 );

      als0 = 4.*PI/beta_0/L2;
      als1 = als0 - pow(als0,2)*b1*LL2;
      als2 = als1 + pow(als0,3)*(pow(b1,2)*(pow(LL2,2) - LL2 -1.) + b2);
      als3 = als2 + pow(als0,4)*(pow(b1,3)*(-pow(LL2,3)+5./2.*pow(LL2,2)+2*LL2-1./2.)-
                              3.*b1*b2*LL2 + b3/2.);

     return als3; 

}
double q_evolution_to_RGI(int Nf,double ainv,double mu_2)
{
        double cmu; // c=cmu
        //mu_2=a2*p2 (adimensional quantity)

        // alphas @ NNLO
        double alm;
        alm=alphas(Nf,ainv,mu_2)/(4*PI);
        ////////////////////////////////
        // N3LO FORMULA
        // Assuming landau gauge
        ///////////////////////////////////
       if(Nf==2){
        cmu = 1. + 2.03448 * alm + 35.9579 * pow(alm,2) + 1199.16 * pow(alm,3);
       }if(Nf==0){
        cmu = 1. + 2.0303 * alm + 42.1268 * pow(alm,2) + 1728.43 * pow(alm,3);
       }if(Nf==4){
        cmu = 1. + 2.4000 * alm + 29.6724 * pow(alm,2) + 719.141 * pow(alm,3);
       }



        return cmu;
}

